#'
#'  Project title:     'Sovereign Risk and Government Change: Elections, Ideology and Experience'
#'  Authors:           Sarah M. Brooks; Raphael Cunha; Layna Mosley
#'  File description:  Analysis of fiscal policy volatility by partisanship and populism
#'  Output:            Table 3; Tables and Figures in Section F of Supplementary Appendix
#'

# my_packages <- c("tidyverse", "countrycode",
#                  "WDI", "rio", "jcolors", "xtable")
# install.packages(my_packages)

library(tidyverse)
library(WDI)
library(countrycode)
library(rio)
library(jcolors)
library(xtable)

# Set directories

DATADIR <- "~/.../Data/Monthly"

# Load EMBIG and CDS spreads dataset

load(file.path(DATADIR, "EMBIG and CDS Dataset.RData"))

# Load WDI data on govt consumption and GDP

wdi <- WDI(country = "all",
           indicator = c("NE.CON.GOVT.KD",  # General govt final consumption expenditure (constant 2010 US$)
                         "NY.GDP.MKTP.KD",  # GDP (constant 2010 US$)
                         "FP.CPI.TOTL.ZG"), # Inflation (annual %)
           start = 1980,
           end = 2015,
           extra = TRUE) 

colnames(wdi)[colnames(wdi) == "NE.CON.GOVT.KD"] <- "gov_consumption"
colnames(wdi)[colnames(wdi) == "NY.GDP.MKTP.KD"] <- "gdp"
colnames(wdi)[colnames(wdi) == "FP.CPI.TOTL.ZG"] <- "inflation"

wdi$iso3c <- countrycode(wdi$iso2c,
                         origin = "iso2c",
                         destination = "iso3c",
                         warn = TRUE)

# Merge partisanship data

dpi <- import(file = file.path(DATADIR, "DPI2015.dta"), format = "dta")
dpi$year <- as.integer(dpi$year)
dpi$iso3c <- countrycode(dpi$ifs,
                         origin = "wb",
                         destination = "iso3c",
                         warn = TRUE)

wdi <- left_join(wdi, dpi[c("iso3c", "year", "execrlc")],
                 by = c("iso3c", "year"),
                 na_matches = "never")

wdi$execrlc <- ifelse(wdi$execrlc == 1, "Right",
                      ifelse(wdi$execrlc == 2, "Center",
                             ifelse(wdi$execrlc == 3, "Left",
                                    ifelse(wdi$execrlc == 0, "No info",
                                                  ifelse(wdi$execrlc == -999, NA, NA)))))

# Add populism data

popul <- read.csv(file.path(DATADIR, "Populist government_v2.1.csv"),
                  stringsAsFactors = FALSE)

popul$iso3c <- countrycode(popul$COWcode,
                           origin = "cown",
                           destination = "iso3c",
                           warn = TRUE)

popul$populist <- popul$patb3
popul <- popul[order(popul$iso3c, popul$year), ]

wdi <- left_join(wdi,
                 popul[c("year", "iso3c", "populist")],
                 by = c("year", "iso3c"),
                 na_matches = "never")

# Interact partisanship and populism

wdi$rlcXpopulist <- ifelse(!is.na(wdi$execrlc) & !is.na(wdi$populist),
                           paste(wdi$execrlc, wdi$populist, sep = " "),
                           NA)
range(unique(popul$year))

# Add info on chief executive from REIGN

reign <- read.csv(file = file.path(DATADIR, "REIGN_2017_01.csv"),
                  stringsAsFactors = FALSE)

reign$iso3c <- countrycode(reign$ccode,
                           origin = "cown",
                           destination = "iso3c",
                           warn = TRUE)

# Leadercode by year

reign <- reign %>%
  group_by(iso3c, year) %>%
  dplyr::filter(month == min(month)) %>%
  dplyr::summarise(leadercode = max(leadercode))

wdi <- left_join(wdi, reign[c("iso3c", "year", "leadercode")],
                 by = c("iso3c", "year"),
                 na_matches = "never")

# Estimate fiscal policy volatility by country

# Log vars

wdi$log_gov_consumption <- log(wdi$gov_consumption)
wdi$log_gdp <- log(wdi$gdp)

# First differences

wdi <- wdi %>%
  group_by(iso3c) %>%
  mutate(d_log_gov_consumption = c(NA, diff(log_gov_consumption)),
         d_log_gdp = c(NA, diff(log_gdp)),
         l_d_log_gov_consumption = dplyr::lag(d_log_gov_consumption, n = 1),
         l_d_log_gdp = dplyr::lag(d_log_gdp, n = 1))

# Limit analysis to period with both partisanship and populism data

wdi <- subset(wdi, year >= 1980) # Same period as populism data

# OLS estimation

wdi$resid <- NA

countries <- unique(wdi$iso3c[!is.na(wdi$iso3c)])

wdi <- subset(wdi, !is.na(iso3c))

for (i in countries){
  
  dat <- wdi[wdi$iso3c == i, ]
  
  # ERROR HANDLING
  possibleError <- tryCatch(
    mod <- lm(d_log_gov_consumption ~ d_log_gdp, data = dat, na.action = na.exclude),
    error = function(e) e)
  
  if(inherits(possibleError, "error")) next
  
  mod <- lm(d_log_gov_consumption ~ d_log_gdp, data = dat, na.action = na.exclude)
  res <- resid(mod)
  wdi$resid[wdi$iso3c == i] <- res
}

# Volatility of fiscal policy by leader

vol_dat <- wdi %>%
  group_by(iso3c, leadercode) %>%
  summarise(fiscal_vol = sd(resid, na.rm = TRUE))

# Summarize partisanship and populism by leader tenure

wdi_leaders <- wdi %>%
  group_by(leadercode) %>%
  summarise(execrlc = first(na.omit(execrlc)),
            populist = first(na.omit(populist)),
            rlcXpopulist = first(na.omit(rlcXpopulist)))

vol_dat <- left_join(vol_dat,
                     wdi_leaders[c("leadercode", "execrlc", "populist", "rlcXpopulist")],
                     by = c("leadercode"),
                     na_matches = "never")

# Policy volatility by govt type

tab_dat <- subset(vol_dat, iso3c %in% wdi$iso3c[wdi$income %in% c("Lower middle income", "Upper middle income")])

tab_countries <- unique(tab_dat$iso3c)

tab1 <- tab_dat %>%
  group_by(execrlc) %>%
  filter(execrlc != "No info") %>%
  summarise(mean_vol = mean(fiscal_vol, na.rm = TRUE),
            sd_vol = sd(fiscal_vol, na.rm = TRUE),
            n = n()) %>%
  mutate(mean_vol = round(mean_vol, 3),
         sd_vol = round(sd_vol, 3))

tab2 <- tab_dat %>%
  group_by(execrlc, populist) %>%
  filter(execrlc != "No info" & !is.na(populist)) %>%
  summarise(mean_vol = mean(fiscal_vol, na.rm = TRUE),
            sd_vol = sd(fiscal_vol, na.rm = TRUE),
            n = n()) %>%
  mutate(mean_vol = round(mean_vol, 3),
         sd_vol = round(sd_vol, 3))

# Print tables

desired_order <- c("Left", "Right", "Center")
tab1 <- tab1[order(match(tab1$execrlc, desired_order)), ]
tab2 <- tab2[order(match(tab2$execrlc, desired_order)), ]

print(xtable(t(tab1)))
print(xtable(t(tab2)))

# Populism by partisanship

my_tab <- table(vol_dat$execrlc, vol_dat$populist)
my_tab

prop.table(my_tab, 1)
xtable(prop.table(my_tab, 1))

# Plot distribution of policy volatility across govt types

plot_dat <- vol_dat

plot_dat$populist <- factor(plot_dat$populist,
                            levels = c("Populist", "Non-populist"))

p <- ggplot(data = plot_dat[-which(is.na(plot_dat$rlcXpopulist) | plot_dat$execrlc == "No info"), ],
            aes(x = fiscal_vol, fill = execrlc)) +
  geom_density(alpha = .5, size = .3) +
  facet_grid(populist ~ execrlc) +
  theme_bw(base_family = "Myriad Pro") +
  theme(panel.grid.minor = element_blank(),
        strip.text = element_text(size = 11, face = "bold"),
        strip.background = element_blank(),
        legend.position = "none") +
  labs(x = "Volatility of fiscal policy", y = "Density") +
  scale_fill_jcolors(palette = "default") +
  NULL
p
